%% Stock Return Extrapolation, Option Prices, and Variance Risk Premium
% By Adem Atmaz
% This file generates Table 5: Effects for the stock: Model vs. evidence
% It also requires 1 function file: olshac.m
% Typically it takes less than 20 seconds to run

clear all;  clc;
format compact
format

% Parameter Values follow from  "Atmaz_Table_4_ParameterValues.m"
alpha=0.50;
theta=0.25;
mu=0.0212;
kappa=0.3567;
zeta=0.0044;
sigma=0.1625;
rho=0;
r=0.0056;

D0=100;
a=2.75;
vec_theta=[0, 0.25, 0.50, 0.75];
N_th=length(vec_theta);


% Initialize vectors
vec_RP=zeros(1,N_th);
vec_vol=zeros(1,N_th);
vec_corr=zeros(1,N_th);
vec_RV=zeros(1,N_th);
vec_VSR=zeros(1,N_th);
vec_BRE=zeros(3,N_th);
vec_BRV=zeros(3,N_th);
for i=1:N_th
    theta=vec_theta(i);    % Varies theta
    % Implied Quantities after fixing theta
    M_DELTA=((kappa+(rho*sigma)*(1+theta))^2)-((2+theta)*(1+theta)*sigma^2);
    M_c=theta/(alpha*(1+theta));
    M_b=-(2+theta)/(kappa+((rho*sigma)*(1+theta))+sqrt(M_DELTA));
    M_sigma_1=(1+(rho*sigma*M_b))/(1-alpha*M_c);
    M_sigma_2=(sqrt(1-rho^2)*sigma*M_b)/(1-alpha*M_c);
    M_Var_con=(M_sigma_1^2)+(M_sigma_2^2);
    M_Cov_con=(rho+(sigma*M_b))/(1-alpha*M_c);
    M_corr_sv=M_Cov_con/sqrt(M_Var_con);
    M_zeta_v=zeta*M_Var_con;
    M_kappa_v=kappa+(sigma*M_Cov_con);
    M_eta_v=(theta/M_b)*(1-(alpha*M_c))*M_Var_con;
    M_sigma_v=sigma*sqrt(M_Var_con);
    
    % long-run mean (steady-state) values of processes
    V0=zeta/kappa; % long run mean of  fundamental  variance
    v0=M_zeta_v/kappa;  % LRM of stock return variance v
    X0=(r+(0.5*M_Var_con*V0))/(1+theta); % LRM of sentiment X
    % Current values of processes
    vt=v0;
    Xt=X0;
    hp_ret=1.20;  % high past return of 20% in past year as in Greenwood
    Xt_change=(1-exp(-alpha))*log(hp_ret);
    vec_X=[X0+Xt_change X0 X0-Xt_change];  %vector of sentiment for high past, SS, low past
    
    % ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    %  ####################### Stock risk premium     ########################
    %  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    vec_RP(i)=vt-(theta*Xt);
    % ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    %  ####################### Stock return volatility ########################
    %  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    vec_vol(i)=sqrt(vt);
    % ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    %  ####################### Correlation ########################
    %  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    vec_corr(i)=M_corr_sv;
    
    % ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    %  #######################   Instantenous  RV       ########################
    %  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    RV_A= M_zeta_v;
    RV_B=1-kappa;
    vec_RV(i)=(RV_A+(RV_B*v0))*(10000/12); %normalized to get monthly values
    
    % ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    %  #######################     Instantenous VSR     ########################
    %  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    VSR_A=M_zeta_v;
    VSR_B=1-M_kappa_v;
    VSR_C=M_eta_v;
    vec_VSR(i)=(VSR_A+(VSR_B*v0)+(VSR_C*X0))*(10000/12); %normalized to get monthly values;
    
    % ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    %  ##################### Bias in return expectation ######################
    %  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    
    for j=1:length(vec_X)
        Xt=vec_X(j);
        vec_BRE(j,i)=theta*Xt;
    end
    % ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    %  ##################### Bias in return volatility #######################
    %  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    for j=1:length(vec_X)
        vt=v0;
        Xt=vec_X(j);
        Subj_var_dt=M_zeta_v+((1-kappa)*vt)+(M_eta_v*Xt);
        Subj_vol_dt=sqrt(Subj_var_dt);
        Obj_var_dt=M_zeta_v+((1-kappa)*vt);
        Obj_vol_dt=sqrt(Obj_var_dt);
        vec_BRV(j,i)=Subj_vol_dt-Obj_vol_dt;
    end
end
vec_VRP=vec_RV-vec_VSR; %Instantenous VRP normalized

disp('Stock risk premium');         round(vec_RP*100,2)
disp('Stock return volatility');    round(vec_vol*100,2)
disp('Correlation');                round(vec_corr,2)
disp('Bias in return expectation'); round(vec_BRE*100,2)
disp('Bias in return volatility');  round(vec_BRV*100,2)





disp('Press any key to continue to Predictive regression ...');
pause;

% ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%  ####################### Predictive regression of D/S ####################
%  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


% Simulation quantities
N_sim=1000;                             % Number of simulations We use 1000
S_end_year=62;                          % Simulation end year. We use 1000
int = 1/252;                            % time interval. Use 1/252 for daily
time_grid = 0:int:S_end_year;           % Time vector from 0 to end year
N_time=length(time_grid);
N_rets=N_time-1;                        % Number of returns in each path
H=1/12;                     % will use monthly returns to run regressions
N_pp=H/int;                 % number of data points for each monthly return

% 1 path vectors
vec_D = zeros(1,N_time); % Allocate output vector
vec_V = zeros(1,N_time); % Allocate output vector
vec_X = zeros(1,N_time); % Allocate output vector
% All paths monthly matrices
N_month=(12*S_end_year)+1;
SM_V = zeros(1,N_month); % Allocate output vector
SM_X = zeros(1,N_month); % Allocate output vector
SM_D = zeros(1,N_month); % Allocate output vector
SM_S = zeros(1,N_month); % Allocate output vector
% Univariate regression of ln(D/S) on 12 month returns
vec_beta_12        =zeros(N_sim,N_th);
vec_R2_12          =zeros(N_sim,N_th);
vec_R2_adj_12      =zeros(N_sim,N_th);
vec_SE_12          =zeros(N_sim,N_th);
vec_SE_hac_12      =zeros(N_sim,N_th);
vec_tstat_12       =zeros(N_sim,N_th);
vec_tstat_hac_12   =zeros(N_sim,N_th);
% Univariate regression of ln(D/S) on 60 month returns
vec_beta_60        =zeros(N_sim,N_th);
vec_R2_60          =zeros(N_sim,N_th);
vec_R2_adj_60      =zeros(N_sim,N_th);
vec_SE_60          =zeros(N_sim,N_th);
vec_SE_hac_60      =zeros(N_sim,N_th);
vec_tstat_60       =zeros(N_sim,N_th);
vec_tstat_hac_60   =zeros(N_sim,N_th);

rng('default');
vec_Z1=randn(N_sim,N_time);
vec_Z2=randn(N_sim,N_time);
tic
for s = 1:N_sim
    s
    for th=1:N_th
        theta=vec_theta(th);
        % Stock price quantities
        M_DELTA=((kappa+(rho*sigma)*(1+theta))^2)-((2+theta)*(1+theta)*sigma^2);
        M_c=theta/(alpha*(1+theta));
        M_b=-(2+theta)/(kappa+((rho*sigma)*(1+theta))+sqrt(M_DELTA));
        M_a=a;
        M_sigma_1=(1+(rho*sigma*M_b))/(1-alpha*M_c);
        M_sigma_2=(sqrt(1-rho^2)*sigma*M_b)/(1-alpha*M_c);
        M_Var_con=(M_sigma_1^2)+(M_sigma_2^2);
        M_Cov_con=(rho+(sigma*M_b))/(1-alpha*M_c);
        % Initalization of each path
        D0=100;
        V0=zeta/kappa; % long run mean of  fundamental  variance
        X0=(r+(0.5*M_Var_con*V0))/(1+theta); %  inital  X is equal to its long run mean
        vec_D(1)=D0;
        vec_V(1)=V0;
        vec_X(1)=X0;
        
        for j = 1:N_rets
            Z1=vec_Z1(s,j);
            Z2=vec_Z2(s,j);
            vec_V(j+1) = max(0, vec_V(j))+((zeta-(kappa*max(0, vec_V(j))))*int)...
                +(sigma*rho*sqrt(max(0, vec_V(j)))*sqrt(int)*Z1)...
                +(sigma*sqrt(1-rho^2)*sqrt(max(0, vec_V(j)))*sqrt(int)*Z2);
            vec_X(j+1) = vec_X(j)+((r+(0.5*M_Var_con*vec_V(j))-(1+theta)*vec_X(j))*alpha*int)...
                +(alpha*M_sigma_1*sqrt(max(0, vec_V(j)))*sqrt(int)*Z1)...
                +(alpha*M_sigma_2*sqrt(max(0, vec_V(j)))*sqrt(int)*Z2);
            vec_D(j+1) = vec_D(j)+((mu*vec_D(j))*int)+(vec_D(j)*sqrt(max(0, vec_V(j)))*sqrt(int)*Z1);
            
        end
        
        % Get end of month values for each process from daily values
        SM_V= vec_V(1:N_pp:end);
        SM_X= vec_X(1:N_pp:end);
        SM_D= vec_D(1:N_pp:end);
        SM_S= SM_D.*exp(M_a+(M_b*SM_V)+(M_c*SM_X)); %stock price
        
        
        % Univariate Regression of 1 year returns on D/S
        % ===============================================
        SM_Ret_12m=(SM_S(12+1:end)./SM_S(1:end-12))-1;
        reg_Y=SM_Ret_12m';
        LT=length(reg_Y);
        X=SM_D(1:end-12)./SM_S(1:end-12);
        reg_X=[ones(N_month-12,1) X'];
        lag=round(0.75*(LT^(1/3)));    %  recommended hac lag choice
        
        [beta, R2, R2adj, X2_NW, X2_HH, X2_R, std_NW, std_HH, std_R, t_NW, t_HH, t_R]=olshac(reg_Y,reg_X,lag,lag);
        
        vec_beta_12(s,th)=beta(2);
        vec_tstat_hac_12(s,th)=t_NW(2);
        vec_R2_12(s,th) = R2;
        vec_R2_adj_12(s,th) = R2adj;
        
        % Univariate Regression of 5 year returns on D/S
        % ===============================================
        
        SM_Ret_60m=(SM_S(60+1:end)./SM_S(1:end-60))-1; %5 year returns
        reg_Y=SM_Ret_60m';
        LT=length(reg_Y);
        X=SM_D(1:end-60)./SM_S(1:end-60);
        reg_X=[ones(N_month-60,1) X'];
        lag=round(0.75*(LT^(1/3)));    %  recommended hac lag choice
        
        [beta, R2, R2adj, X2_NW, X2_HH, X2_R, std_NW, std_HH, std_R, t_NW, t_HH, t_R]=olshac(reg_Y,reg_X,lag,lag);
        
        vec_beta_60(s,th)=beta(2);
        vec_tstat_hac_60(s,th)=t_NW(2);
        vec_R2_60(s,th) = R2;
        vec_R2_adj_60(s,th) = R2adj;
        
        
    end   % end theta
    
end    % end sim


clc;
toc

%  REPORT UNIVARIATE REGRESSION for 12 months returns
beta_12m=mean(vec_beta_12);
tstat_hac_12m=mean(vec_tstat_hac_12);
R2_12m=mean(vec_R2_12);
disp('Beta for 12m returns'); round(beta_12m ,2)
disp('t-stat(hac) for 12m returns');  round (tstat_hac_12m,2)
disp('R2 for 12m returns in %'); round(R2_12m*100,2)


%  REPORT UNIVARIATE REGRESSION for 60 months returns
beta_60m=mean(vec_beta_60);
tstat_hac_60m=mean(vec_tstat_hac_60);
R2_60m=mean(vec_R2_60);
disp('Beta for 60m returns'); round(beta_60m,2)
disp('t-stat(hac) for 60m returns'); round(tstat_hac_60m,2)
disp('R2 for 60m returns in %'); round(R2_60m*100,2)


clear vec_Z1 vec_Z2 Z1 Z2 time_grid ;
clear vec_D vec_S vec_V vec_X;
%% THE END













